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Abstract: 

We study the behaviour of the normal derivative of eigenf unctions of the Helmholtz equation 
inside billiards with Dirichlet boundary condition. These boundary functions are of particular 
importance because they uniquely determine the eigenfunctions inside the billiard and also 
other physical quantities of interest. Therefore they form a reduced representation of the 
quantum system, analogous to the Poincare section of the classical system. For the normal 
derivatives we introduce an equivalent to the standard Green function and derive an integral 
equation on the boundary. Based on this integral equation we compute the first two terms of 
the mean asymptotic behaviour of the boundary functions for large energies. The first term is 
universal and independent of the shape of the billiard. The second one is proportional to the 
curvature of the boundary. The asymptotic behaviour is compared with numerical results for 
the stadium billiard, different limagon billiards and the circle billiard, and good agreement is 
found. Furthermore we derive an asymptotic completeness relation for the boundary functions. 
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1 Introduction 



The study of eigenf unctions of quantum systems, in particular with chaotic classical dynamics, 
has attracted a lot of attention. A prominent class of examples is provided by billiard systems, 
which classically are given by the free motion of a particle inside some domain with elastic 
reflections at the boundary The corresponding quantum system is described by the Helmholtz 
equation inside a compact domain flcR 2 (in units h = 1 = 2m), 

Atfj n (x) + k 2 n ij) n (x) = , x G fi , (1) 

with (for example) Dirichlet boundary conditions 

if> n (x) = , x G dQ , (2) 

where the normalized eigenfunctions ip n (x) are in L 2 (Q). A detailed knowledge about the 
behaviour of eigenvalues and the structure of eigenstates is relevant for applications, for example 
in microwave cavities or mesoscopic systems (see e.g. [|IJ and references therein). 

A particular nice feature of the classical dynamics in Euclidean billiards is the existence of a 
global Poincare section based on the boundary dQ. This Poincare section facilitates the study 
of the classical dynamics considerably, because the dynamics is reduced to an area-preserving 
map on the two-dimensional compact surface of section. It appears therefore desirable to look 
for a similar reduced representation of the quantum mechanical problem. One representation is 
provided by the boundary integral method. The boundary integral method transforms the two- 
dimensional Helmholtz equation (|l|) for the eigenfunctions with Dirichlet boundary condition 
to a one-dimensional integral equation on the boundary dfl. This method involves the normal 
derivative of the eigenfunctions, hereafter called the boundary functions, 

d 

u n {s) := -^-ipn(x)\ x=x{s) = (n(s),Vip n (x(s))) , (3) 

where x(s) is a point on the boundary dQ, parameterized by the arclength s, and n(s) denotes 
the outer normal unit vector to dQ at x(s). The integral equation for the boundary function 
is of the form 

u = -H x {k)u , (4) 

where Hi(k) is an integral operator depending on the parameter k (for the explicit form of its 
integral kernel see eq. fl25|) below). This equation has solutions only for a discrete set of values 
of the parameter k which, when real valued, give the eigenvalues of the Helmholtz equation. 
The associated solution u is the normal derivative of the corresponding eigenfunction which 
can be obtained from u via an integral formula using the free Green function on the plane. 

This reduction to the boundary is very useful for the numerical computation of eigenvalues 
and eigenfunctions. The boundary function allows for a direct expression of the corresponding 
eigenfunction, its normalization [FJ], momentum distribution ||, autocorrelation function || 
and other quantities of interest. Furthermore the boundary functions are the basis to define a 
Husimi representation of the eigenstates over the classical Poincare section (see e.g. 0, |6|) and 
therefore provide a direct connection with the classical Poincare map. This is in particular useful 
in situations where one is interested in fine structures of the eigenstates and their relation to the 
classical dynamics, like in the field of quantum chaos. From this it is clear that the boundary 
functions deserve a study in their own right. A profound knowledge of their properties can then 
be used to obtain a description for the above mentioned quantities. 
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There are two main aspects concerning the boundary functions. First one can consider 
these functions as a possible set of basis vectors that span a kind of natural space for the 
reduced quantum system. Here the question of orthogonality and completeness of the boundary 
functions arises. Second, as the boundary functions contain all information on the quantum 
system inside the given domain, it is very interesting to see how the properties of the eigenstates 
of the Helmholtz equation ([I]) are reflected in the boundary functions. The last point is our 
main interest in this paper where we study the mean semiclassical behaviour of the sequence 
of boundary functions in terms of a spectral average. 

A classical example of a spectral quantity is the spectral staircase function (integrated level 
density) 

N(k) := #{n E IN | k n < k} , (5) 
whose asymptotic behaviour for k — > oo is given by the Weyl formula (see [0]) 

N{k) = ±k>-±k + o{h) . (6) 

Here A denotes the area of the billiard and £ the length of the boundary dQ. It is a well known 
observation that the first two terms of (||) usually describe the mean behaviour of N(k) very 
well, even down to the ground state (see for example fig. 2 in f8j). 

In the same way one can consider the sum over a sequence of the normalized eigenfunctions 
ip n up to some given energy k 2 , 

®(k,x) := \Mx)\ 2 • (7) 
For billiards with C°°-boundary || Theorem 17.5.10] implies 

where d(x) is the distance of the point x e Q to the boundary dfl, and J\{z) denotes the Bessel 
function. The remainder R(k, x) satisfies the estimate \R(k,x)\ < Ck for large k. The second 
term in (^j) describes the influence of the boundary. This result is, as the Weyl formula for 
the mean behaviour of the spectral staircase function, completely independent of the classical 
chaoticity of the underlying system. So it applies equally well to integrable, mixed and chaotic 
systems. To illustrate the behaviour of ^f(k,x), let us consider a member of the family of 
limagon billiards introduced by Robnik 0, |TTJ, whose boundary is given in polar coordinates 
by p(v 9 ) — 1 + £ cos<y9 where e G [0, 1] is the family parameter. In fig. [1] we show a three- 
dimensional plot of ty(k,x)/k 2 for the desymmetrized limagon billiard with parameter e = 0.3 
using the first 100 eigenfunctions of odd symmetry (i.e. with Dirichlet boundary conditions on 
the symmetry axis). For this parameter value the classical billiard has a mixed phase space 
(see fig. 1 in f!2]| ). The agreement of the asymptotic behaviour given by (|[) with ty(k,x) is 



very good (see [ 13 1 for further examples). 

In analogy a similar behaviour of a sum over the sequence of boundary functions u n (s) is 
expected in the limit k — > oo, i.e. 

<k n 
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Figure 1: Three-dimensional plot of ^(k,x)/k 2 for the desymmetrized limagon 
billiard for e = 0.3 and k chosen such that N(k) = 100. To the right a section of 
ty(k,x)/k 2 at x\ = 0.5 for the case N(k) = 1000 is shown and compared with the 
asymptotic result (dashed curve). The horizontal dotted line shows the leading 
term 1/(4tt). The inset shows a magnification. 



In |14| it has been conjectured that the asymptotic behaviour of a similar sum, J2k <k \ u n(s)\ 2 , 
is Dik A + D 2 n(s)k 3 , where k(s) is the curvature of the boundary at the point s G dQ. In this 
work we derive for the sum (§) the asymptotic behaviour c\k 2 + C2ff(s)/c including the constants. 
The first constant turns out to be c\ = 4- which is consistent with the leading term of the Weyl 
formula (|B|). Multiplying (^) by ^(n(s),x(s)) and integrating over the billiard boundary gives 
N(k) for the left hand side, because of the following normalization relation for u n (s) (for 



alternative derivations and more general boundary conditions see [ IE , IB ) 



1 

2 

an 

For the right hand side one uses 



(n(s),x(s))\u n {s)\ 2 ds = k 2 n . (10) 



i j(n(s),x(s))ds = A (11) 



no 



to obtain the leading term of the Weyl formula The next to leading term can be determined 
by using 

J K(s)(n(s),x(s))ds = £ , (12) 



an 



which for billiards with C°°-boundary follows from 4-t(s) = —k(s)ti(s) and partial integration 
(here t(s) is the unit tangent vector at s G dfl). 

In the following we derive the full asymptotic series for the sum in eq. @ and compute the 
first two terms explicitly. To this end in sec. El an integral equation on the boundary is derived 
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which is then used in sec. [3] to obtain the mean behaviour of the boundary functions. Some of 
the more technical details are given in the appendices. In sec. |] we provide a numerical test of 
the asymptotic behaviour at finite energies. We conclude with a short summary. 



2 An integral equation on the boundary 

Our main tool for the analysis of the semiclassical behaviour of the boundary functions is a 
boundary Green function which we define as 

£la^^(«K0O- ( 13 ) 
new n 

This is the analogue of the standard Green function inside the billiard, where we have replaced 
the eigenf unctions by their normal derivatives on the boundary. Our aim in this section is to 
derive an integral equation for fllcf), which can be solved recursively, and therefore leads to an 



expansion of (13) in terms of known functions. The method is very similar to the one used 
by Balian and Bloch [FT], and our expansion on the boundary is the exact analogue of their 
multiple reflection expansion. 

The energy-dependent Green function for the Dirichlet Laplacian on Q is given by 

G(£,*,y) = £^fM (14) 



new 



for (x, y) £ Vt x Q, where k 2 and ip n are the eigenvalues and normalized eigenfunctions of A, 
respectively, and E is a complex parameter. This Green function is the unique solution of the 
equation 

(A x + E)G(E,x,y) = 5(x-y) (15) 

in Q which satisfies Dirichlet boundary conditions. To obtain a similar relation on the boundary 
we define 

d d 

g(k, s, s') :=——G(k 2 , x, y)\ x = x (s), y = y (s') 

= Z) fc 2_ fc 2 M "(sK(5') , 

new n 

where k = \J~E~ 6 C, Imfc > 0, is the branch of the square root of E which has positive 
imaginary part. 

We will now derive an integral equation for g(k, s, s'). Let 

Go(E, x, y) := JL J __L_ e ^> d 2 £ (17) 

be a free Green function which satisfies (A x + E)Gq(E, x, y) = 5(x — y) on R 2 . This function 
depends holomorphically on E in the cut plane C\R + . Additionally we introduce the auxiliary 
functions f{x,s') := -£^G(E, x, y)\ y = y ( s ') and f {x,s') := -^G (E,x,y)\ y=y ( s ,y Then 

(A x + E)(f(x,s')-f o (x,s')) = (18) 
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for x G Q \ dQ and therefore we can represent this difference as a single layer potential 

f(x, s') - f (x, s') = SL^x, s') := J G (E, x, y(s))fi(s, a') ds . (19) 

The density /jl(s, s') is determined by the boundary condition 

-Mx(s"), s') = J G (E, x(s"), y(s)Ms, s') ds , (20) 
an 



and it is a standard result from potential theory that this equation is solvable fL9| . We can 
determine /x(s, s') from the jump relations for a single layer potential JO 



d± x SL,(x(s),s') = J d niB G (E,x(s),y(s"))ti(s",s / ) ds" ± ^(s, s') , (21) 



mi 



where denotes interior (+) or exterior (— ) limits of the normal derivative, i.e., for a function 
4>{x) they are defined as d^<p(x(s)) := \im E ^ (n(s),'V x (f)(x(s) =pen(s))) where n(s) denotes 
the outer normal unit vector to the boundary at x(s). From equation (|T9| ) and the jump 
relations ( F2T|) we then obtain 

g(k, s, s') = d+J(x(s), s') - d-J(x(s), s') = fi(s, s') . (22) 

Applying now to equation (JHfy and using Q2T]) and (f2^) leads to the desired integral equation 
on the boundary 

g(k, s, s') = g {k, s, s') - J hi(k, s, s")g{k, s", s') ds" , (23) 

an 

where g(k,s,s') is given by (|llj) and 

g (k, s, s') = 2d nx d ny G (k 2 , x{s), y(s')) , (24) 
hx{k, s, s') = 2d nx G (k 2 , x{s),y{s')) . (25) 

Equation (f23|) is a Fredholm integral equation of second kind and can be solved by iteration. 
If we write G,Go,Hi for the operators with integral kernels g,go and hi, respectively, the 
integral equation becomes 

G = G — HiG , 

which can be solved for G as 

G = (1 + HiY x G = i-lTH^Go , (26) 

nSINo 

and the series converges if Hi is small enough. Going back to the kernels of the operators this 
gives the expansion 

oo „ 

g(k, s, s') = g (k, s, s') + ^T(-l) n / h n (k, s, s")g (k, s", s') ds" (27) 

n=1 an 
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where for n > 2 we have 

h n (k, s, s") = J J ■ ■ ■ J hi(k,s, st)hi(k, s 1 ,s 2 )--- h(k, s n _ 1; s") ds 1 ds 2 • • • ds n _i . (28) 
an an an 

Let us now discuss the convergence of the series (|26|). By the integral representation (|64| ) 
for hi(k,s,s') (derived in appendix A) we see that for s ^ s' a positive imaginary part Imk 
gives an exponential damping factor, therefore we expect that there exists a 7 > such that 
for Im k > 7 the operator norm of H 1 satisfies 

H-ffill < 1 . (29) 

Then, for Im k > 7, the series ( p6|) would converge to a bounded operator. On the other hand, 
the boundary integral method tells us that the operator Hi has an eigenvalue — 1 if k 2 is an 
eigenvalue of the Helmholtz equation (p]). At these values the series ( p6|) clearly diverges, as it 
should, since g has a pole there. 

It is often useful to consider spectral functions different from g(k,s,s'), and one way of 
introducing them such that the series ( pT|) can still be applied is based on the formula 



1 2(7; + 

lim Im - — 1 T, , „ = S(x -k) + 5(x + k) . (30) 

7^0+ 7T (x + 17) 2 - k 2 

Let p(k) and a(k) be functions which are holomorphic in the strip — e < Im k < 7 + e for some 
e > 0, real valued for real arguments and a(k) even. Furthermore assume that their product 
a(k)p(k) decays faster than 1/k 2 for Re A; — > ±00. Then we define 

00+17 

g p (k, s, s') := Im / Izpik — z)a(z)g(z, s, s') dz . (31) 

7T J 

— oo+i7 



Under the assumption ( p9|) for Im k > 7 we can insert the expansion ([27]) and obtain a series 
which converges to the kernel of a bounded operator. On the other hand, since g p (k, s, s') is 



holomorphic in {k E C | Imk > 0}, the integral (|3T| ) does not depend on 7. By taking the limit 



7 — > + and using (|30|) and the definition (|Tj|) we get the representation 

9 p (k, s, s') = Y,lp( k - k n) + P(k + k n )]a{k n )u n {s)u* n {s') . (32) 
new 

So the role of the function p is to select a spectral window in the summation, whereas a(k) acts 
as a weight function. 

Now p{k) and a(/c) can be adapted to the particular question one is interested in. In the 
next section we want to study the mean behaviour of the boundary functions, and therefore we 
will choose 

P(k) = P(t)e itk dt (33) 

where p is even and has compact support in a sufficiently small neighborhood of 0. This choice 
obviously fulfils the requirements needed. For the weight function a(k) we will choose 

«(*) = (^) 

with a > 7. This also fulfils the requirements and satisfies a(k) = 1/k 2 + 0(l/k ) for k — > 00, 
which gives the correct normalization factor for the u n in view of (|T0D. 
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3 Mean behaviour of boundary functions 

In this section we want to study the mean behaviour of the boundary functions u n for large 
energies k 2 ,. We will do this by choosing a suitable test function p in fl3lD , namely we will 
assume that there is an e > such that 

supp/5 C [—£,£] and p(t) = 1 for t G [— e/2, e/2] , (35) 

and furthermore that e is smaller than the shortest periodic orbit of the classical billiard flow. 
From fl2"?D , fl28|) and (|3l|) we obtain the expansion 

g p (k,s,s')= J29 P n(k,s,s') (36) 

nelo 

with 

oo+i7 

(-I)" 4 " 1 f f 

gp(k,s,s'):= Im / 2zp{k — z)a(z) / h n (z, s, s )go(z, s , s ) ds dz (37) 

— oo+i7 dQ 

for n > 1, and for n = 

oo+i7 

g£(k, s, s') := Im / 2zp(k — z)a(z)g Q (z, s, s') dz . (38) 

—00+17 

If p is furthermore positive, the sum (^) can be interpreted as defining a mean value of 
boundary functions weighted with the factor a(k n ), where the mean is taken over a spectral 
window around k defined by p. In the following we will assume that the billiard boundary is 
smooth, or can be obtained as the desymmetrization of a smooth billiard. We will show in 
appendix |C| that under the conditions (|35]) and with a{k) ~ 1/k 2 for k — > 00 

g^k,s,s)=0(k 1 - n ) , (39) 

and therefore the sum (|36|) provides an asymptotic expansion for large k of the mean boundary 
functions. 

The explicit computation of the first terms in the expansion ( pB]) is given in appendix [B|. 
Here we choose a to be of the form (0). For s — s' the first term is given by 

gg(k,s,s) = ^ + 0(k-°°) , (40) 



and for s ~ s' one has 

k r 2 

9o(k, s, s') — — 



k\s — s 

The second term is for s ~ s' given by 



-Mkls-s'D + Ois-s') 



+ 0(AT°°). (41) 



9i{k, s, s') = cos(2A;| S - a '|) + 0(1/*) , (42) 

Z7T 
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where k(s) denotes the curvature of the boundary at s. 
Integrating g p (k, s, s') one obtains 



A- 

(N*p)(k,s,s>) = JgP(k',s,s')dk' (43) 



where * denotes convolution and N(k, s, s') is defined as 

m ,s,sr.= ^ U -^^. (44) 

It is well known from the theory of spectral asymptotics that 

N(k, s, s') = (N*p) (k, s, s')(l + 0(l/k)) , (45) 

see 



and therefore from an asymptotic expansion of we immediately obtain the leading 
semiclassical behaviour of ([44]) . 

Inserting (|36|) and using (|40|), (|42l) and ( [39l) the asymptotic behaviour of (|43T) becomes 



(N*p)(k,s,s) = ^k 2 -^k + 0(\nk) , (46) 

which is the main result of this section. 
Thus the leading term of (|45|) reads 

N(k,s,s) = ^k 2 + 0(k) . (47) 

Assuming in addition that the set of p £ [—1,1], such that (p,s) belongs to a periodic orbit 
of the billiard map, has measure zero, we expect that the two-term asymptotics ( [46]) holds for 
N(k,s,s) as well, but with an error term o(k). However, to prove this requires to adapt the 
more sophisticated methods from [I8| or j7|, where similar statements are proven for N(k). 



Note that, as discussed in the introduction, the relation between N(k) and N(k, s, s) is given 
by 

N(k) = l -J (n(s), x(s))N(k, s, s) ds , (48) 
an 

and therefore we recover (Q) from ((IB] ) using fllUp, ( ]TT| ) and (|T2]), up to the error term. 

Often one studies billiards Q with discrete symmetries and restricts the study to the cor- 
responding symmetry subclasses of the eigenfunctions. For example for a system which is 
symmetric with respect to reflection at the Xi-axis the eigenfunctions can be classified as either 
odd, fulfilling i/}(xi,X2) = —ip(xi,—X2), or even, where ip(xi,X2) = ip(x\, — x 2 ). Consequently 
eigenfunctions with odd symmetry satisfy Dirichlet boundary conditions on the symmetry axis 
and even eigenfunctions obey Neumann boundary conditions on the symmetry axis. Of course, 
such symmetries of Q induce symmetries of dQ. We restrict ourselves to the case of a reflection 



symmetry at a point s £ <9f2. In this case the boundary Green function (|T3|) has to be modified 
to 

g ± (k, s, s') = g(k, s + s, s + s') ± g(k, s + s, s - s') (49) 
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for s, s' near s with even (+) or odd (— ) symmetry. Using this and ( PD we obtain for s close 
to s 

I ± - f Ji(2fc|s - sol) 



k 



g^(k,s,s) = — 
Ztt 



k\s 



s 



0{l/k)) 



and 

Therefore we get 

{N ± *p)(k,s,s) 



K[S) 

' 2tt 



[1 ± cos(2fc|s - sol)] (1 + 0(l/k)) . 



47T 



1 ± 



1- J (2k\s - s \) 
\s — s | 2 fc 2 



MA) 

' 2tt 



1 ± 



sin(2A;|s — sq\ 
2k\s — sq| 



(50) 
(51) 

+ 0(\nk) , (52) 



for s close to a fixpoint So of the symmetry. 

Our results show that the mean behaviour of the normalized boundary functions is very 
similar to the mean behaviour of eigenfunctions. The crucial difference between the two se- 
quences of functions {ip n }n£K, {u n }n£K is that the eigenfunctions live on a two-dimensional 
space whereas the boundary functions live on a one-dimensional space. Since both u n and ip n 
oscillate roughly with the same de Broglie wave length 2ix fk w this leads to an overcompleteness 
of the set {w n }n£iN- This statement can be made more explicit by observing that (fill) implies 



gP(k,s,s') = -5(s-s') + 0(l/k) . 

71 



More precisely, this means that for every if e C°°(dQ) 



0(l/k) 



holds with coefficients 



7T 

2^1 



(53) 



(54) 



(55) 



an 



This follows from (|72|) by the method of stationary phase (see appendix |D|). Since p is a rapidly 
decreasing function, this means that the boundary functions with spectral parameter k n in an 
interval of fixed width around k form a complete set in the limit k —>■ oo. The number of these 
states grows like k, in contrast to the number of all states up to energy k 2 , which, according 
to the Weyl formula, grows like k 2 . Therefore this result gives a quantitative measure of the 
overcompleteness of the set {u n } ne -^. 



4 Numerical results for integrable, mixed and chaotic 
systems 

In this section we would like to test how well the asymptotic expansion (fiBD describes the mean 
behaviour of the boundary function at finite energies. To illustrate the energy dependence we 
plot in fig. |2] N(k, s, s) for s = 3.2 for the desymmetrized limagon billiard with e = 0.3 (see the 
introduction). The dashed line beneath the staircase function (full curve) is the asymptotic 
result (f4Tf). As for the spectral staircase function excellent agreement, even down to the ground 
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Figure 2: Plot of N(k, s, s) at s — 3.2 for the desymmetrized limagon billiard with 
e = 0.3. The dashed curve is the asymptotic result fl46|) and the dotted curve is just 
the leading term (pFTj). 



state, is observed. Notice, that for such a good agreement it is crucial to include the curvature 
term. This is illustrated by the dotted line in fig. || which is a plot of the leading term k 2 / (47r). 

Now we turn to the s-dependence of N(k, s, s), where for better comparison we have divided 
by k 2 



a(k, 



-y 

k 2 ^ 



\u, 



k 2 



pN(k,s,s) 



(56) 



First we consider the stadium billiard which is given by two semi-circles joined by two parallel 
straight lines. The stadium billiard is proven to be strongly chaotic, i.e. it is ergodic, mixing 
and a i^-system pD| , |2"T|| . The height of the desymmetrized billiard is chosen to be 1, and a 
denotes the length of the upper horizontal line, for which we have a = 1.8 in the following. The 
boundary is parameterized starting with s = at the corner of the quarter circle, ranging to 
s = 7r/2 at the place where straight line and quarter circle join tangentially until s = 7r/2 + etat 
the next corner. Fig. ^ shows a(k, s) using the first 1000 boundary functions. The asymptotic 
result (see sec. 0) reads (for even and odd symmetry) 



1 ± 



E 



1 - J (2k\s - Si \ 



corners i 



M5J 

2nk 



1 ± 



E 



sin(2fc|s — Sj|) 
2k\s — sA 



corners i 



(57) 



where the contributions near the corners are taken into account. 

At s = tt/2 the mean behaviour shows a jump caused by the discontinuity in the curvature 
of the boundary (transition from the circular part to the straight line). Although our derivation 
is not valid at this point we observe excellent agreement of a(k, s) with a~(k, s). Moreover, the 
behaviour near the two corners is clearly visible and very well described by W~(k, s). 
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Figure 3: Plot of the spectral average a(k, s) for the fully desymmetrized stadium 
billiard with a = 1.8 and Dirichlet boundary conditions everywhere. The energy k 2 
is chosen such that the first 1000 boundary functions are taken into account. The 
horizontal dotted line is the leading term, l/(4n), and the dashed line corresponds 
to the asymptotic formula a~(k,s), eq. (|57|) . The inset shows a magnification. 
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Figure 4: Same as in the previous figure but for the desymmetrized limagon billiard 
with e = 0.3. For this parameter the classical phase space is mixed. 
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Figure 5: Same as in the previous figure but for the desymmetrized limagon billiard 
with e — 0.7. Here the influence of the curvature contribution is clearly visible. 
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Figure 6: Plot of a(k, s) for the desymmetrized circle billiard with radius one and 
Neumann boundary condition compared with the corresponding asymptotic result 
a + (k,s) shown as dashed line. 
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The next example, shown in fig. |] for the limagon billiard with e = 0.3, illustrates that the 
mean behaviour is independent of the classical dynamical properties as for this parameter the 
considered billiard is a mixed system (i.e. regular and irregular regions in phase space coexist). 
The asymptotic formula also works well in cases where the curvature changes more strongly as 
for example for the limagon billiard with e = 0.7, see fig. [5[ In this case the classical dynamics 
appears to be ergodic (though it is not the case since there exist some very small islands of 
stability ||). 



Finally, we would like to consider the circle billiard as an example of an integrable system. 
The boundary functions of a circle billiard with radius R and the corresponding eigenvalues are 
given by 

3 lm J-l-% 7„2 _ ( Jim 



Ulm{s) = ~7^ e R and kL = {W ' (58) 

respectively, where ji m denotes the m-th positive zero of the Jz-Bessel function. Therefore, 

km(s)| 2 1 



N(k) . (59) 



Turning to the desymmetrized circle billiard no such simple expression exists. Fig. || shows a 
plot of a(k, s) for the desymmetrized circle billiard with Neumann boundary condition on the 
symmetry axis. Again we find very good agreement with the corresponding asymptotic result 
o+(k,s). 



5 Summary 

We have studied the semiclassical behaviour of the normal derivative of eigenfunctions of Eu- 
clidean billiards. These boundary functions form a reduced representation of the quantum 
system, analogous to the Poincare section of the classical billiard. Thus they are of special in- 
terest in order to understand the properties of the eigenfunctions of billiards and their relation 
to the classical system, as investigated in the field of quantum chaos. 

We have introduced an analogue of the standard energy Green function for the boundary 
functions, and one of our main results is the derivation of an integral equation for this bound- 
ary Green function. The resulting integral equation is the principal tool for the semiclassical 
analysis and leads to an expansion of the boundary Green function which is the analogue of the 
classical multiple reflection expansion. The resulting expansion could now be used to derive a 
semiclassical representation of the boundary Green function in terms of orbits of the classical 
billiard map, but here we were only interested in the so-called length-zero contribution which 
determines the mean behaviour of the normal derivatives for large energies. 

For the mean behaviour of the sequence of the boundary functions we have derived a two- 
term asymptotics for large energies, where the first term is universal, i.e. completely independent 
of the given billiard, and the second term is proportional to the curvature of the boundary. 
These theoretical results fit very well with the numerical computations for different Euclidean 
billiards. 

Furthermore we have obtained a completeness relation for the boundary functions with 
momenta k n near k for large k, which shows that the whole set of normal derivatives is highly 
over-complete, since a fraction of order k is asymptotically sufficient to span the reduced state 
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space over the boundary. 
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Appendix 

A Integral representations 

For the computation of the two leading terms in the expansion (|36"D, eqs. fl40|) and (0), we 
need suitable integral representations of hi(k,s,s') and g (k,s,s'). We start by deriving a 
representation of hi(k,s,s'). Namely, from (|i~7|) and (|25|) we get 

hi(k, s, s') = J^~y J ^§ d 2 £ , (60) 



where E = k 2 . Inserting the relation 

1 



|2 



is / e< E -\^ dt , (61) 



o 



E-\£\ 

with e = sgnlmi? and for | Im£?| > 0, allows us to solve the ^-integral and we arrive at 

oo 

hi(k, s, s>) = -| J M°)> S *M) e H^¥^l dt ( (62) 



where we have used the abbreviation 5x(s, s') = x(s) — y(s'). This expression has the disad- 
vantage that the factor (n(s), 8x(s, s')) tends to zero for s — > s', so we use 

1 d • ifa( s ,a')i 2 ie \so>(3, s ')\ 2 

;e 1£ ^^- = tt^^^ (63) 



\5x(s,s')\ 2 dt At 2 

and partial integration (the boundary term at t = vanishes in the weak sense as a function 
of s) to obtain 

oo 

V^y, = _^ (n(»),fa(». y)) f eH ^ + ^ d( . (64) 
vr |5aj(s, s')r y 



Notice that for s ~ s' 

{n{a) ' SxM) -U^+0{s-^ t (65) 



|&r(s,s')| 2 2 

where k(s) denotes the curvature of the boundary dfl at x(s) (with the sign convention that 
it is positive for a circle). 
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A similar representation can be derived for g (k,s,s f ). Inserting now ( |5"T| ) in ( |2~3| ) allows 
again to compute the ^-integral and we obtain 



go(k,s,s') 



2tt 



(n(s), 5x(s, s'))(n(s'), 5x(s, s')) 



4r 3 

. (n(s),n(s')} 
16 2r 2 



(66) 



The integral is again well defined as an oscillatory integral, but it will be useful below to have 
a suitable regularisation at hand when we want to change the order of integration. To this end, 
let Go,({k 2 , x, y) be the integral kernel of the complex power of the free resolvent, (— A + k 2 )^^, 
( G C, and let go,((k, s, s') := 2d nx d riy Go^{k 2 , x(s), y(s')). Then a calculation analogous to the 

one leading to (|66|), using (E — |^| 2 ) - '' = e r ^ £ / r C- 1 et 1£ ( E -\£\ ) r d r instead of (|61|), shows that 



£ C+l e -if(C-l) 



27rr(C) 



(n(s), 5x(s, s'))(n(s'), 5x(s, s')) 



4r 3 

. (n(s),n(s')) 



(67) 



r (-l e ie[±(6 X (s,s>)) 2 +k 2 r] dr 



and we have of course <?o,i(^> s , s ') — 9o(k, s, s'). 



B The computation of and g p x 

In this appendix we will use the integral representations derived in appendix [X| to compute the 
first and second term in the asymptotic expansion (|36|). 

We start with the first term, (38). Since the integrand is holomorphic in the upper half- 
plane, we can take the limit 7 — > 0, and using (Ell) we obtain 



9 P o(k,s,s') = j^yj +P(*+ l*l)MI*l)^ d 2 £ . (68) 

R 2 

Since p is rapidly decreasing, and |£| is positive, the term p(k + |£|) is smaller than any 
negative power of k for large k and can therefore be neglected. Using the representation (pI3| ) 
and introducing polar coordinates in the ^-integral leads to 

2n 00 00 

2k A 



g P o(k,s,s') = —J J j p(t)a(kr)(n(s),e( V ))(n(s'),e( V )) 

-00 

el k[t(l-r)+r(5 X (s,s')M^)] ^ dtdrd^ + 0(k~°°) , 

where e((p) denotes the unit vector in direction if. The r, t-integrals can now be evaluated and 
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give 

oo oo 

h I I me imi ~ r)+r{5x(s ' S ' )Mv))] a(rk)r 3 dtdr (70) 

-oo 

oo oo 

= A J J p(t+ (5x(s, s r ), e(<p)))e- iktr a((r + l)k)(r + 1) drdt e^i'^Mv)) 

— oo —1 

= e^p(t + (6x{s, s'),e(tp)))a{{r + l)k){l + r)\ r=t=0 e ik ^ s > s '^)) + 



for | (5a; (s, s'), e(ip))\ < e/2 by (p5|) . In the step from the second to the third line we have used 
the general relation 

oo oo 

7^ J I e- ikrt f(r,t)drdt = e-^f(r,t)\ r=0 , =0 . (71) 

— oo — oo 

The error 0(k~°°) is due to the cut-off of the r-integral at r = —1. This finally leads to 

3 

2tt 2 " 





Z7T 

g p (k,s,s') = / (n( S ),e(v,))(n( S '),e(^))e ife ^')'^ dip + 0(Ar«) , (72) 



for s close to s'. At s = s' we get 



a, 5) = + 0(k~°° ) , (73) 

Z7T 



and for s ~ s' a Taylor expansion of the exponent gives 



2tt 



g£(k,a,s' 



a(k)k 3 f 2 ifc |W| s in0 



2tt 2 



cos e 

(74) 
27T k|s — s'l 

For the determination of the second term, g^(k, s, s'), it is useful to use the regularized 
expression ( |67|) with ( in a range where all integrals converge, and finally make an analytic 
continuation to ( = 1. So we use the representations (|64| ) and ([67]) and insert them in fl3"7|). 
The ^-dependence in the expressions for g ,c an d ^i is simple and we can perform the resulting 
^-integral 

oo+i7 

— / 2zp(k — z)a(z)hi(z, s, s")go r(z, s" , s') dz . (75) 

IT J 

—00+17 

Since the integrand is holomorphic in the upper half-plane, we can perform the limit 7 — > 0. 
If we furthermore change variables r — > r/k and t — > in ( |66D and ([641) , and interchange the 
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order of integration, the ^-integral boils down to 

00 00 
i J 2z 3 p(k-z)a(z)e it -^ z2 dz = ^-e ik{t+r) J 2(k - zfa{k - z)p(z)e it ^ L * e -a<t+r> dz 
—00 —00 



(76) 



Here we need that a satisfies a kind of symbol estimate. Collecting the remaining integrals 
leads to 



9i X (k,s,s') = Im 



k 4 ^a(k) 



IT" 



an 



b(k,s,s', s",r,()p(2(t + r))(l + 0(l/k)) 



e ik[^(S X (s'',s')f + ±(8 X (s,s'')) 2 + (r+t)] ds "(irdt 



(77) 



with 



b{K,s,s,s ,r,y e \ 5x (s,s")\ 2 



' (n(s"),6x(s" , s'))(n(s'), 5x(s", s')) 

1/C 



+ 



4^3 

(n(s"),n(s')) 



2r 2 



(78) 



This looks quite complicated and to understand the properties of this expression better let 
us discuss the stationary points of the phase function ^(5x(s" ', s')) 2 + -^(5x(s, s")) 2 + (r + t) 
with respect to s",t,r. The derivatives with respect to t and r give the conditions 



-^(^ S ")) 2 + i = o , -^(M*V)) 2 + i = o 

respectively, and hence 2t = \Sx(s, s")\ and 2r = \5x(s", s')\. The ^''-derivative leads to 

±(t(s"),5x(s",s')) - ±(t(s»),6x(s,s»)) = 
which yields together with the previous conditions on t and r for s" 

(t(0,«i(?>)> = (t(s"),5x^7")} . 



(79) 



(80) 



(81) 



But this is just the condition that there exists a trajectory in the billiard Q starting at the point 
s at the boundary, which is then elastically reflected at s" and ending in s'. This is of course 
what one expects by analogy with similar expressions, namely that g p (k, s, s') is semiclassically 
given by a sum over all classical orbits from s to s', each contributing an amplitude depending 
on the stability of the orbit and an oscillating factor with frequency proportional to the length 
of the orbit. The n-th term in the expansion of g p (k, s, s') contains exactly the orbits with 
n reflections on the boundary. For the determination of the contribution of these orbits in 
leading order one can simplify the formulas for g and hi considerably by using their asymptotic 
expansions for large arguments which can be easily derived by the method of stationary phase 
from the integral representations in appendix A. But if s is close to s', there is one very short 
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orbit with s" between s and s' whose length tends to zero for s — > s'. The contribution of 
this orbit determines the mean behaviour of g p (k, s, s) and therefore it is commonly called 
the length-zero contribution. For the computation of this contribution the above mentioned 
asymptotic formulas for go and h\ cannot be used because they are not valid in these regions, 



hence we must work with the full representation as in (|77|) . 

The ^''-integral can be solved by the method of stationary phase, and since p is supported 
in a small neighborhood of there is only one stationary point for s close to s', and we find 
that g[^{k, s, s') is for s ~ s' equal to the imaginary part of 

' ' / \ e I * l ^T +c * +r)J p(2(t + r)) drdt (1 + 0(l/fc)) . 

(82) 



27r 5 / 2 r(C) J J r 5 /2-C^t + r 





Introducing now the coordinates v — t + r and w = t — r gives 

e ik[L W^T +it+r)] p(2(t + r)) drdt 



OO OO . 

1 f t V' 2 ; ,r(£^)! 






oo 



2^ [v-W I } V + W Xi\ dw e^^^p(2v) dv 
J J (v - wfl 1 "^ 







(83) 



OO 

= B (~, + J v- 3 / 2+ <e ik ^ + ^p(2v) dv 

x ivrp(O) 2 e-^i-^iJ^C^Is - + 0(1/A:)] , 

where -B(it, t>) denotes the beta function and H„(x) the Hankel function of the first kind. After 
collecting all terms we can perform the analytic continuation to £ = 1 and finally obtain 



gi (k, s, s') = -p(0)a(k)k 2 ^- cos(k\s - s'\)[l + 0(l/k)} (84) 



for s ~ s'. 



C Estimating g p n 

In this appendix we derive the estimate ( ff9"D on s, s). For g?(k, s, s) one obtains a similar 
expression as (|77|), 



CO 



gUk,s,s) =Im [k n+l J J j b(s 1 s\r,r)a(r/k)p(2r + 2j2t^ 

o k™ an« 1=1 (85) 
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with 



71+1 



b(s, s',r) = Yl 



n si 



«0> 



1^(4 



(n(s' n ), 5x(s' n , s)){n(s), 6x(s' n , s)) 



4r 3 



— l- 



: (n(s' n ),n(s)) 



2r 2 



(86) 



where s' = (s[, ■ ■ ■ , sJJ, t = (ti, ■ ■ ■ ,t n ) and we use the convention s' = s n+ i = s. As in the 



discussion after (|77|) one sees that the main contributions to (|85| ) come from orbit segments, 
starting in s and returning to s after n reflections at the boundary. The total length of these 
orbits is 2r + 2 £\ £j and since this expression appears as the argument of p, the only orbit 
contributing to the integral d85|) is the one with s[ = s' 2 = ■ ■ ■ = s' n = s and t« = r = 0, thanks 
to the small support of p. So we can approximate (8x(s' i _ 1 , s^)) 2 ~ (s-_ x — s-) 2 for all z in the 
exponent of the integrand of (|85|) , and in the prefactor b we can set s- = s for all i. Now a 
substitution s' — > s'/k, t t/k and r — > r/A; makes the integrand independent of up to a 
factor of /c _2n , which together with the former prefactor A; n+1 gives an overall factor k~ n+1 . A 
more careful analysis shows that the resulting integrals diverge at r = 0, so we should use the 
regularization which we already applied to the computation of gi(k,s,s'). But this does not 
change the final result 

g n (k,s,s)=0(k 1 - n ) . (87) 



D The completeness relation 

In this appendix we derive the completeness relation (Q). By the results of appendix |C| we get 

J y{s')g p {k, s, s') ds' = J ^(s')g p (k, s, s') ds' + J ^s'^k, s, s') ds' + 0(1/*) (88) 

an an an 

if (p G C°°(dQ) and p satisfies the conditions (|55|). The first term on the right hand side can 
be easily computed by using (|72|), the method of stationary phase, and by observing that due 
to the cutoff introduced by p the only stationary points come from s' = s. The result is 



J <p(s')g p (k } s, s') ds' = 2 - V (s) + 0(1 A) 



(89) 



<-»> 



where we have assumed that a(z) = 1/z 2 + 0(l/z 4 ), see (Bl). The computation of the second 



term in (88) is similar but more complicated. Using now (|77| ) and solving the resulting s' and 
^''-integrals with the method of stationary phase leads to 



^(^g^k^s.s') ds' 



an 



Im 



Ck 3 ^K(s)cp(s)e-^ (<:+1) 




t l/2 
-5/2-< 



p(2(t + r))e ifc( * +r) drdt 



(90) 



;i + 0(l/fc)) , 



o o 



where we have collected all factors not depending on k and s in the real constant C . Introducing 
new coordinates v = t + r,w = t — r allows to solve the integrals, and after setting ( = 1 we 
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arrive at 



J V {s')g[{k, s, s') ds' = C"^ tp(a) + 0(1/ k 2 ) = 0(1/ k) (91) 
an 

for k — > oo, with another constant C" . Therefore the result (0) is established. 
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